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Abstract. Wc study the computational complexity certification of inexact gradient augmented 
Lagrangian methods for solving convex optimization problems v^ith complicated constraints. We solve 
the augmented Lagrangian dual problem that arises from the relaxation of complicating constraints 
with gradient and fast gradient methods based on inexact first order information. Moreover, since 
the exact solution of the augmented Lagrangian primal problem is hard to compute in practice, we 
solve this problem up to some given inner accuracy. We derive relations between the inner and the 
outer accuracy of the primal and dual problems and we give a full convergence rate analysis for both 
gradient and fast gradient algorithms. We provide estimates on the primal and dual suboptimality 
and on primal feasibility violation of the generated approximate primal and dual solutions. Our 
analysis relies on the Lipschitz property of the dual function and on inexact dual gradients. We also 
discuss implementation aspects of the proposed algorithms on constrained model predictive control 
problems for embedded linear systems. 
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1. Introduction. Embedded control systems has been widely used in many ap- 
plications and its usage in industrial plants has increased concurrently. The concept 
behind embedded control is to design a control scheme that can be implemented on 
autonomous electronic hardware, e.g a programmable logic controller [24], a micro- 
controller circuit board [25l [20] or field-programmable gate arrays [8]. One of the 
most successful advanced control schemes implemented in industry is model predic- 
tive control (MPC) and this is due to its ability to handle complex systems with hard 
input and state constraints. MPC requires the solution of an optimal control problem 
at every sampling instant at which new state information becomes available. In the 
recent decades there has been a growing focus on developing faster MPC schemes, 
improving the computational efhciency |18| and providing worst case computational 
complexity certificates for the applied solution methods [9l [HI [TU [20] , making these 
schemes feasible for implementation on hardware with limited computational power. 

For fast embedded systems [71 [HI [10] the sampling times are very short, such that 
any iterative optimization algorithm must offer tight bounds on the total number of 
iterations which have to be performed in order to provide a desired optimal controller. 
Even if second order methods (e.g. interior point methods) can offer fast rates of 
convergence in practice, the worst case complexity bounds are high [2]. Further, 
these methods have complex iterations, involving inversion of matrices, which are 
usually difficult to implement on embedded systems, where the units demand simple 
computations. Therefore, first order methods are more suitable in these situations (9] 
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When the projection on the primal feasible set is hard to compute, e.g. for 
constrained MPC problems, an alternative to primal gradient methods is to use the 
Lagrangian relaxation to handle the complicated constraints and then to apply dual 
gradient schemes. The computational complexity certification of gradient-based meth- 
ods for solving the (augmented) Lagrangian dual of a primal convex problem is studied 
e.g. in [U m HU [121 [m [171 [TH [22] . In |1] the authors present a general framework for 
gradient methods with inexact oracle, i.e. only approximate information is available 
for the values of the function and of its gradient, and give convergence rate analysis. 
The authors also apply their approach to gradient augmented Lagrangian methods 
and provided estimates only for dual suboptimality. In |22| an augmented Lagrangian 
algorithm is analyzed using the theory of monotone operators. For this algorithm 
the author proves asymptotic convergence under general conditions and local linear 
convergence under second order optimality conditions. In |171ll6j a dual fast gradient 
method is proposed for solving quadratic programs with linear inequality constraints 
and estimates on primal suboptimality and infeasibility of the primal solution are pro- 
vided. In }llj the authors analyze the iteration complexity of an inexact augmented 
Lagrangian method where the approximate solutions of the inner problems are ob- 
tained by using a fast gradient scheme, while the dual variables are updated by using 
an inexact dual gradient method. The authors also provides upper bounds on the 
total number of iterations which have to be performed by the algorithm for obtaining 
a primal suboptimal solution. In [12] a dual method based on fast gradient schemes 
and smoothing techniques of the ordinary Lagrangian is presented. Using an averag- 
ing scheme the authors are able to recover a primal suboptimal solution and provide 
estimates on both dual and primal suboptimality and also on primal infeasibility. 

Despite widespread use of the dual gradient methods for solving Lagrangian dual 
problems, there are some aspects of these methods that have not been fully studied. 
In particular, the previous work has several limitations. First, the focus is mainly 
on the convergence analysis of the dual variables. Second, only the dual gradient 
method is usually analyzed and using exact information. Third, there is no full con- 
vergence rate analysis (i.e. no estimates in terms of dual and primal suboptimality and 
primal feasibility violation) for both dual gradient and fast gradient schemes, while 
using inexact dual information. Therefore, in this paper we focus on solving convex 
optimization problems (possibly nonsmooth) approximately by using an augmented 
Lagrangian approach and inexact dual gradient and fast gradient methods. We show 
how approximate primal solutions can be generated based on averaging for general 
convex problems and we give a full convergence rate analysis for both algorithms that 
leads to error estimates on the amount of constraint violation and the cost of primal 
and dual solutions. Since we allow one to solve the inner problems approximately, 
our dual gradient schemes have to use inexact information. 

Contribution. The contributions of this paper include the following: 

1. We propose and analyze dual gradient algorithms producing approximate pri- 
mal feasible and optimal solutions. Our analysis is based on the augmented 
Lagrangian framework which leads to the dual function having Lipschitz con- 
tinuous gradient, even if the primal objective function is not strongly convex. 

2. Since exact solutions of the inner problems are usually hard to compute, we 
solve these problems only up to a certain inner accuracy ein. We analyze 
several stopping criteria which can be used in order to find such a solution 
and point out their advantages. 

3. For solving outer problem we propose two inexact dual gradient algorithms: 
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- an inexact dual gradient algorithm, with complexity 0(l/eout) iterations, 
which allows us to find an ^out-optimal solution of the original problem by 
solving the inner problems with an accuracy £;„ of order 0{eont)- 

- an inexact dual fast gradient algorithm, with complexity ©(-^/l/eout) itera- 
tions, provided that the inner problems are solved with accuracy ejn of order 

C'(£out-\/^out)- 

4. For both methods we show how to generate approximate primal solutions and 
provide estimates on dual and primal suboptimality and primal infeasibility. 

5. To certify the complexity of the proposed methods, we apply the algorithms 
on linear embedded MFC problems with state and input constraints. 

Paper outline. The paper is organized as follows. In Section (TJ motivated by 
embedded MFC, we introduce the augmented Lagrangian framework for solving con- 
strained convex problems. In Section[2]we discuss different stopping criteria for finding 
a suboptimal solution of the inner problems and provide estimates on the complexity 
of finding such a solution. In Section [3] we propose an inexact dual gradient and fast 
gradient algorithm for solving the outer problem. For both algorithms we provide 
bounds on the dual and primal suboptimality and also on the primal infeasibility. In 
Section |4] we specialize our general results to constrained linear MFC problems and 
we obtain tight bounds on the number of inner and outer iterations. We also provide 
extensive numerical tests to prove the efficiency of the proposed algorithms. 

Notation and terminology. We work in the space R" composed by column vec- 
tors. For x,t/ G R", {x,y) :— x^y — ^-'i^/i and := (X^iLi -^D^^^ denote the 
standard Euclidean inner product and norm, respectively. We use the same notation 
(•,•) and II • II for spaces of different dimension. We denote by conejoi, i € /} the 
cone generated from vectors {a^, i G /}. We also denote by i?p :— max^^yg^ \\z — y\\ 
the diameter, int(Z) the interior and hd{Z) the boundary of a convex, compact set 
Z. By dist(j/, Z) we denote the Euclidean distance from a point y to the set Z and 
by hz{y) := sup^g^y'^^ the support function of the set Z. For any point z G Z we 
denote by Mz{z) := {s \ {s,z~z)<Oyz€ Z} the normal cone of Z at z. For a real 
number x, [x\ denotes the largest integer number which is less than or equal to x, 
while ":=" means "equal by definition". 

1.1. A motivating example: Linear MPC problems with state-input 
constraints. We consider a discrete time linear system given by the dynamics: 

Xk+i = A^Xk + BuUk, 

where Xk G M""' represents the state and Uk G M"" represents the input of the system. 
We also assume hard state and input constraints: 

XkGX C M"- , UkGU C M"" Vfc > 0. 

Now, we can define the linear MPC problem over the prediction horizon of length N , 
for a given initial state follows [231: 

s.t. x^+i = A^Xi + BuUi, xq = X, 
Xi e X, ui € U Vz, e Xf, 

where both the stage cost £ and the terminal cost £f are convex functions (possibly 
nonsmooth). Note that in our formulation we do not require strongly convex costs. 
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Further, the terminal set X{ is chosen so that stability of the closed-loop system is 
guaranteed. We assume the sets X, U and to be compact, convex and simple (by 
simple we understand that the projection on these sets can be done easily, e.g. boxes). 
Furthermore, we introduce the notation z := \^x^ ■ ■ ■ u"^ ■ ■ ■u^_^ , Z := 

riila^ ^ X Xf X Y\a=i U ^i^d /(^) '■— J2iLo^ ^i^iT^i) + ^{{^n)- We can also write 
compactly the linear dynamics Xi+i ~ A^Xi + BuUi for all i = 0, • • • , — 1 and 
xo = a; as Az = b{x) (see [23l[26j for details). Note that b{x) € R^"- depends hnearly 
on X, i.e. b{x) := [{A^x)'^ 0^ • • • O"'"] . In these settings, for linear MPC we need to 
solve, for a given initial state x, the primal convex optimization problem: 

(P(a;)) mm{f{z) \ Az ^b{x), z e Z}, 

z 

where / is a convex function (possibly nonsmooth) and A is a matrix of appropriate 
dimension. Moreover, the set Z is simple as long as X, X{ and U are simple sets. In the 




following sections we discuss how we can efficiently solve optimization problem (P(a;) I 
approximately with dual gradient methods based on inexact first order information 
and we provide tight estimates for the total number of iterations which has to be 
performed in order to obtain a suboptimal solution in terms of primal suboptimality 
and infeasibility. 

1.2. Augmented Lagrangian framework. Motivated by MPC problems, we 
are interested in solving convex optimization problems of the form: 

(p) r 

where / is convex function (possibly nonsmooth), A G ]|J'"X" is a full row-rank matrix 
and Z is a simple set (i.e. the projection on this set is computationally cheap), compact 
and convex. We will denote problem ((P|) as the primal problem and / as the primal 
objective function. 

A common approach for solving problem ((P|) consists of applying interior point 
methods, which usually perform much lower number of iterations in practice than 
those predicted by the theoretical worst case complexity analysis On the other 
hand, for first order methods the number of iterations predicted by the worst case com- 
plexity analysis is close to the actual number of iterations performed by the method 
[14) . This is crucial in the context of fast embedded systems. First order methods ap- 
plied directly to problem (|P]) imply projection on the feasible set {z \ z G Z, Az = b}. 
Note that even if Z is a simple set, the projection on the feasible set is hard due to the 
complicating constraints Az ~ b. An efficient alternative is to move the complicating 
constraints into the cost via Lagrange multipliers and solve the dual problem approx- 
imately by using first order methods and then recover a primal suboptimal solution 
for ((P|) . This is the approach that we follow in this paper: we derive inexact dual 
gradient methods that allow us to generate approximate primal solutions for which 
we provide estimates for the violation of the constraints and upper and lower bounds 
on the corresponding primal objective function value of ((P|) . 

First let us define the dual function: 

(1.2) d{X) :=imnC{z,X), 

where C{z, A) := f{z) + (A, Az — b) represents the partial Lagrangian with respect to 
the constraints Az = b and A the associated Lagrange multipliers. Now, we can write 
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the corresponding dual problem as follows: 
(D) max d(X). 

We assume that Slater's constraint qualification holds, so that problems ([P|) and ((D|) 
have the same optimal value. We also denote by z* an optimal solution of ((P]) and 
by A* the corresponding multiplier (i.e. an optimal solution of ((D|) ). 

In general, the dual function d is not differentiable [1] and therefore any subgradi- 
ent method for solving has a slow convergence rate [13] . We will see in the sequel 
how we can avoid this drawback by means of augmented Lagrangian framework. We 
define the augmented Lagrangian function [S] : 

(1.3) £p{z,X) ■.^f{z) + {X,Az-b) + ^\\Az-bf, 

where p > represents a penalty parameter. The augmented dual problem, called 
also the outer problem, is defined as: 

(Dp) max dp (A), 

where dp{X) := min£p(z. A) and we denote by 2*(A) an optimal solution of the inner 
problem min^g^ £p(z. A) for a given A. It is well-known [1] [11] that the optimal 



value and the set of optimal solutions of the dual problems (|D|) and (Dpi coincide. 
Furthermore, the function dp is concave and differentiable and its gradient is |15) : 

Vdp(A) ■.^Az*{X)-b. 

Moreover, the gradient mapping V(ip(-) is Lipschitz continuous with a Lipschitz con- 
stant U] given by: 

In conclusion, we want to solve within an accuracy Eout the equivalent smooth outer 



problem (Dpi by using first order methods with inexact gradients (e.g. dual gradient 
or fast gradient algorithms) and then recover an approximate primal solution. In other 
words, the goal of this paper is to generate a primal-dual pair (z. A), with z G Z, for 
which we can ensure bounds on dual suboptimality, primal infeasibility and primal 
suboptimality of order Eout, i-c: 

(1.4) /* - dp(A) < O (eout) , \\Az- b\\ < O (Cout) and |/(z) - /*| < O (cout) • 

2. Complexity estimates of solving the inner problems. As we have seen 
in the previous section, in order to compute the gradient Vrfp we have to find, for a 
given A, an optimal solution of the inner convex problem: 

(2.1) z*(A) e argmin/:p(z. A). 

From the optimality conditions |21| , we know that a point z* (A) is an optimal solution 
of (P?T]) if and only if: 



(2.2) 



(V/:p(z*(A),A),z-z*(A)) >0 VzeZ. 
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An equivalent way to characterize an optimal solution z*{X) of (j2.1l) can be given in 
terms of the following inclusion: 

(2.3) Oe V£p(z*(A),A)+AAz(^*(A)). 

Since an exact minimizer of the inner problem (j2.ip is usually hard to compute, 
we arc interested in finding an approximate solution of this problem instead of its 
optimal one. Therefore, we have to consider an inner accuracy which measures 
the suboptimality of such an approximate solution for (|2.ip : 

^(A) w argmin f/(z) + (A, Az - b) + ^\\Az - b\\^] . 

Since there exist several ways to characterize an £in-optimal solution [H [TTJ [22] , we 
will further discuss different stopping criteria which can be used in order to find such 
a solution. A well-known stopping criterion, which measures the distance to optimal 
value of (|2.ip . is given by: 

(2.4) z{X)eZ, Cp{z{X),X)~Cp{z*{X),X)<el. 

The main advantage of using ()2.4[) as a stopping criterion for finding z(A) consists 
of the fact that in the literature [13] there exist explicit bounds on the number of 
iterations which has to be performed by some well-known first or second order methods 
to ensure the Ein-optimality. 

Another stopping criterion, which measures the distance of z{X) to the set of 
optimal solution Z*{X) of p.ip is given by: 

(2.5) z{X) e Z, dist (z(A), Z*{X)) < O (£;„) • 

It is known that this distance can be bounded by an easily computable quantity when 
the objective function satisfies the so-called gradient error bound [1]. Thus, we can use 
this bound to define stopping rules in iterative algorithms for solving the optimization 
problem. Note that gradient error bound assumption is a generalization of the more 
restrictive notion of strong convexity. 

As a direct consequence of the optimality condition (|2.2p , one can use the following 
stopping criterion: 

(2.6) z{X)eZ, {S/Cp{z{X),X),z-z{X))>-0{ein) "^z e Z. 
Note that (|2.6p can be formulated using the support function as: 

hz i-WCp (z(A), A)) + (V£,(z (A), A) , z(A)) < O (£;„) . 

When the set Z has a specific structure (e.g. a ball defined by some norm), tight upper 
bounds on the support function can be computed explicitly and thus the stopping 
criterion can be efficiently verified. 

Based on optimality conditions ()2.3|) . the following stopping criterion can also be 
used in order to characterize an £i,i-optimal solution z(A) of the inner problem ()2.1[) : 

(2.7) z(A) e Z, dist(0, V£p(z(A), A) +Nz{z{X))) < 0(£i„). 

The main advantage of using this criterion is given by the fact that the distance in 
p.7|) can be computed efficiently for sets Z having a certain structure. Note that 
()2.7|) can be verified by solving the following projection problem: 

(2.8) s* e argmin (z(A), A) + 

seAfziziX)) 
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To see how (|2.8p can be solved efficiently we are first interested in finding an explicit 
characterization of the normal cone J\fz{z(X)), when the set Z has a certain structure. 

Lemma 2.1. Assume that the set Z is a general polyhedral set, i.e. Z := 
{z e R" I Cz < c}, with C G RP^" and c e Rp. Then, problem can be recast as 

the following quadratic optimization problem: 

(2.9) mm\\VCp{z{X),\) + C^^if, 

where matrix C contains the rows of C corresponding to the active constraints in 
Cz(X) < c. In particular, if Z is a box in M", then problem (|2.9p becomes separable 
and it can be solved explicitly in 0{p) operations, where p represents the number of 
active constraints in Cz(A) < c. 

Proof. Let us recall that if z(A) e mt{Z), then we have A/z (z(A)) = {0} and there- 
fore the distance dist (0, VCp (z(A), A) + Afz (z(A))) will be equal to ||V£p (z(A), A) ||. 
In the case of ^(A) € hd{Z) there exists an index set I(z(A)) C {1, • • • ,p} such 
that Ciz{X) = Ci for all i £ I(z(A)), where Ci and Ci represent the i-th row and 
i-th element of C and c, respectively. Using now Theorem 6.46 in [5T] we have 
Afz (^(A)) = cone {Cf , i Gl (z(A))}. Introducing the notation C for the matrix whose 
rows are Ci for alH G I (z(A)), we can write (|2.8p as (|2.9p . Note that, in problem (|2.9p . 
the dimension of the variable ^ is p ~ |I (z(A))| (i.e. the number of active constraints) 
which usually is smaller in comparison with the dimension n of problem p.8|) . 

Now, if we assume that Z is a box in R", then problem (|2.9p can be written in 
the following equivalent form: 

min Ifi'^CC'^fi + VCl (^(A), A) C'^ 

Since for box constraints we have CC'^ = Ip, i.e. the identity matrix, the previous 
optimization problem can be decomposed into p scalar projection problems onto the 
nonnegative orthant and thus in order to compute the optimal solution of (|2.9p we 
only have to perform 0{p) comparisons. □ 

The next lemma establishes some relations between stopping criteria (|2.4p - p.7p . 

Lemma 2.2. The conditions (|2.4p . (|2.5p . (|2.6p and (|2.7p satisfy the following: 

(i) Let V£p be Lipschitz continuous with a Lipschitz constant Lp. Then 

(EH) dm), ([23]) dm) (12:61) . 

(ii) //, in addition, Lp is strongly convex with a convexity parameter Up > 0, then 

(1^ ^ (EH) ^ (1^ . 

Proof (i) (123]) => (1231): In Section 3 of [4] the authors show that if ([231) holds, 
then (1231) also holds with ©(ein) = ef^ + y2L^i?pein. 
(|231) =^ ([231): We can write: 

{VCp{z{\),\),z-z{\)) 

= (V£p( z(A), A )-V£p(z*(A), A) , z-z{\)) + {V Cp{z*{\), A) , z-z*{\) + z*{\) - z(A)) 
> -(Lpi?p + ||V/:p(z*(A),A)||)ein. 

Since Z is compact and V£p(-,A) is continuous, then V£p(-,A) is bounded. Hence, 
our statement follows from the last inequality. 
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(1231) ^ ^M- The condition ^1} can be written as (V£p( z(A), A ), z - z(A)) > 
— £in (e, z — z{X)) for all e such that ||e|| < 1 and z £ Z. If Z is bounded, then the 
last inequality implies (V£p( z(A), A ), z — z{X)) > —RpEm- 

(ii) (|2.7[) (|2.4p : Since £p(z,A) is strongly convex and z*{\) is its minimizer over 
Z ^ we have: 

> C, (z(A), A) - (z*(A), A) > ^||z(A) - z*(A)|p. 

From the convexity of Cp{z, A) we can write: 

£p(z(A), A) - /:p (z*(A), A) < {VCp (z(A), A) , z(A) - z*(A)) 
< (V£p (z(A),A) + s*,z(A)-z*(A)) < || V/Ip (z(A), A) + ||z(A) - z*(A)|| 



-{Cp (z-(A),A)-£p(z*(A),A)) 



1/2 



which implies (|2.4p . 

(|2.4p => (|2.5p : Taking into account that £p(-, A) is strongly convex we have ^ ||z*(A) — 
z(A)|| < £p(z(A),A)-/:p(z*(A),A) < ef„. This leads to |jz(A)-z*(A)|| < (2/crp)i/2ei„. 
The lemma is proved. □ 

The next theorem provides estimates on the number of iterations that are required 
by fast gradient schemes to obtain an ein approximate solution for inner problem p.ip . 

Theorem 2.3. Assume that Junction Cp{-,X) has Lipschitz continuous gra- 
dient w.r.t. variable z, with a Lipschitz constant Lp and a fast gradient scheme 
is applied for finding an ein approximate solution z(A) of such that stopping 

criterion ()2.4p holds, i.e. £p(z(A),A) — £p(z*(A),A) < ef^. Then, the complexity of 

finding z(A) is O (^/^^ iterations. If, in addition Cp{-,-) is strongly convex with a 

convexity parameter Op > 0, then z(A) can be computed in at mo 
iterations by using a fast gradient scheme. 



Note that if the function / is nonsmooth, we have a complexity O ( ^ ^ ] itera- 



tions with a projected subgradient method or an improved O by using smoothing 

techniques [15| . provided that / has a certain structure. 

3. Complexity estimates of solving the outer problem using approx- 
imate dual gradients. In this section we solve the augmented Lagrangian dual 



problem ( Dp ) approximately by using dual gradient and fast gradient methods with 
inexact information and derive computational complexity certificates for these meth- 
ods. Since we solve the inner problem inexactly, we have to use inexact gradients and 
approximate values of the augmented dual function dp defined in terms of z(A), i.e. 
we introduce the following pair: 

dp{X) := Cp (z(A), A) and Vdp(A) := Az{X) - b. 

The next theorem, which is similar to the results in [4], provides bounds on the dual 
function when the inner problem (|2.ip is solved approximately. For completeness we 
give the proof. 

Theorem 3.1. // z(A) is computed such that the stopping criterion (|2.6p is 
satisfied, i.e. z(A) G Z and min^g^ (V£p(z(A), A), z — z(A)) > — (l + yj2LpRp) Si-a, 
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then the following inequalities hold: 

dp{X) + {Vdp{X),fl- A) - - Af - (l + y^Rp) £in < dpi^i) 

(3.1) _ , _ ^ ^ ^ ^ 

<dp{X) + {Vdp{X),fi-X) VA,/ieE'". 



Proof. For simplicity we introduce the notation Cz := 1 + y/2LpRp. The right- 
hand side inequality follows directly from the definitions of dp and dp . We only have 
to prove the left-hand side inequality. Following the derivations from Section 3.3 in 
[4] we have: 

dpifi) > min f /(z(A)) + (V/(z(A)), z-z(A)) + {fi,Az - b) + ^\\Az - bf} 

> mm\f{z{X))~(A'^X + pA^{Az{X)-b),z~z{X)) + {fi, Az - b) + ^\\Az - b\\'^] 
+ min (V£p (z(A), A) , z - z(A)) , 

where we use the convexity of / and the properties of minimum in the first and the 
second inequality, respectively. Using now the assumptions of the theorem and the 
definition of dp we obtain: 

dpin) 

>mm\f{z{X)y{A'^X + pA'^{Az{X)~b),z~z{X)) + {fi,Az-b) + ^\\Az-b\\'^}-Czein 
min{(A(z-z(A)),A*-A) + pA(z-z(A))f}+dp(A)-f(vJp(A),A*-A)-Cz£i„ 



zez 



>min {{A{z~-z{X)),^i-X) + ^\\A(z-z{X))f]+dp{X) + {^dp{X),^l-X)-Czeu 



By taking into account that min^gK" ^IICIP + 9^£. = — ^llfflP and using the definition 
of Cz, we obtain from the last expression the right-hand side inequality. □ 

Note that the first inequality helps us construct a quadratic model which bounds 
from below the function dp, when the exact values of the dual function and its gradi- 
ents arc unknown. The second inequality can be viewed as an approximation of the 
concavity condition on dp. The two models, the linear and the quadratic one, use only 
approximate function values and approximate gradients evaluated at certain points. 
A more general framework for inexact gradient methods can be found in [4] . 

It can be easily proved that under the assumptions of Theorem 13. 1[ the following 
relation holds between the true and approximate gradient: 




II Vdp(A) - Vdp(A)|| < W2id 1 + v/2ipi?p £i„ VA e 



3.1. Inexact dual gradient method. In this section we provide the conver- 
gence rate analysis of an inexact dual gradient ascent method. Let {aj}j>o be a 

sequence of positive numbers. We denote by Sk := '^'j^o'^i- ^^^^ section we 
consider an inexact gradient method for updating the dual variables: 



(IDGM) Afc+i := Xk + afcVJp(Afe), 
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where offe e [i ^, ^] C (0, +00) is a given step size and L > L^. We recall that the 
inexact gradient is given by: 

VJp(Afc) Azk-b, 

where '■= z(Xk). The following theorem provides estimates on dual suboptimality 
for the scheme pPGMp . 

Theorem 3.2. Suppose that the conditions of Theorem \3.1\ are satisfied. Let 
{-^*;}fc>o ^ sequence generated by (jIDGMp and {Xk}k>o be an average sequence 
derived from {Afe}^>Q as Afe := S^^ X]j=o '^j^j+i- Then, for any k > 1, the following 
estimate for the dual suboptimality holds: 



(3.2) 



2{k + l) 



where we define Rd '■= \\Xa — A*||. 

Proof. Let rj := — A*||. By using pPGMp and the estimates p.ip we have: 

r-f+i =r] + 2 (Aj+i - A^, A^+i - A*) - ||Aj+i - A^-f 

- 2a, (VJp(A,), A* - A,) - (1 - a,L^) || A,+i - 



2a i 



La 



(VJp(A,),A,.+i-A,-)-^||A,.+i-A,-! 



< 



(3.3) 



< 



r| + 2a, [dp{X,) - dp{X*)] + 2a, K(A,+i) - d,{X,)] 
+ 2ajCze\n - (1 - ajid)||Aj+i - A^jp 
r] - 2a, [dp{X*) - dp{X,+i)] + 2a,Czen,. 



Here the last inequality follows from a, ^\L_ ^ ,L^^\. Summing up the last inequality 
from J = to fc and taking into account that dp (A*) = /*, we obtain: 

fe 

^ 2a,\f* - dp(A,+i)] < i?^ + 2^feCzein. 
Now, by the concavity of dp and the definition of Afe, this inequality implies: 



Sk 



f*-dpiXk) 



Note that /* - dp{Xk) > and Sk > L^^{k + 1). The last inequality together with 
the definition of Cz imply p.2p . □ 

Next, we show how we can compute an approximate solution of the primal prob- 
lem ((P|) . For this approximate solution we estimate the feasibility violation and the 
bound on the suboptimality for (|P|). Let us consider the following average sequence: 



(3.4) 



Since Zj S Z for all j > and Z is convex, then z^ G Z . From the iteration of 
Algorithm pPGMp and p.4[) . by induction, we have: 



(3.5) 



Afe+i = Ao + Sk{Azk - b). 
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The following two theorems provide bounds on the primal infeasibility and the primal 
suboptimality, respectively. 

Theorem 3.3. Under assumptions of Theorem lS.Sl the sequence Zk generated by 
p.4|) satisfies the following upper bound on the infeasibility for primal problem (jP|) ; 



(3.6) \\Azk-b\\<i^{k,ein), 



where. (Ke..) :^ ^ + J^^^^^'^S^' 



Proof. From (|3.3p we have: 

||A,+i-A*i|2 < \\Xj-X*\\^-2a,[dp{X*)-dp{X,+i)]+2a,Czein < \\Xj-X*f + 2ajCzSin. 

Here the last inequality follows from the fact that dp{X*) — dp{Xj+i) > and aj > 0. 
By induction, it follows from the above inequality that: 

(3.7) ||Afe+i - A*|p < ||Ao - A*||2 + 2SkCzein. 



Now, from p.Sp we have: 

||Afc+i - A*||2 = ||Ao - A* + SkiAzk - b)\\^ > [||Ao - A*|| - Sk\\Az - b\\] 
= ||Ao - X*f - 2Sk\\Xo - X*\\\\Azk - b\\ + SlWAzk - bf. 

Substituting this inequality into (|3.7p we obtain: 

Sl\\Azk^b\\^ - 2Sk\\Xo-X*\\\\Azk-b\\ < 25feCz£i„. 
The last inequality implies: 

_ ^ Rd + [Rl + 2SkCzeinV/^ < 2i?d , f2C7e 



Sk Sk V Sk 

Note that Sk> L^^{k + 1). This inequality implies □ 

Theorem 3.4. Under the assumptions of Theorem \3.3[ the primal suboptimality 
can be characterized by the following lower and upper bounds: 

II All + ^^^(fc,£in)] ^^(fc,ei„) < fizk) -f*< + (l + V^Rp) £in- 



2{k + l) 

Proof. Let us first prove the left-hand side inequality. Since /* = tip (A*), by using 
the definition of Cp{zk, A*) and the Cauchy-Schwartz inequality we get: 

/* - dp{X*) < Cpizk,X*) = f{zk) + {\\Azk -b) + ^\\Azk - bf 

< f{zk) + \\X*\\\\Azk - b\\ + ^\\Azk -bf< f{zk) + \\X*\Hk,e,^) + ^v{k,e,^f. 

Here, the last inequality follows from Theorem 13.31 In order to prove the right-hand 
side inequality we first use the convexity of Cp and the assumptions of Theorem 13.11 

Cp{zj, Xj) < Cp {z*{Xj), Xj) - (\'Cp{zj,Xj),z*{Xj) - zj) < dp{Xj) + CzEin- 
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Previous inequality together with the definition of Cp and dp{Xj) < f* lead to: 

f{z,) + {X„Az, -b) + P-\\A-z, -h\\^- r < Cze,r.. 
Using the iteration of Algorithm pDGM[) and aj < p = L^^ we obtain: 



< 



^(l|A,f-||A,- 



Multiplying this inequality with aj and then summing up these inequalities from j = 
to k we get: 

k 

Now, using the definition of Zk and the convexity of / we can deduce that: 



-r< 



2Sk 



Cz£m- 



Finally, by taking into account that 5*^ > L{k + 1), we get from the last estimate the 
right-hand side inequality. □ 

Let us fix the outer accuracy Eout ■ We want to find the number of outer iterations 
/sout and a relation between Sout and Ein such that after this number of iterations of 
Algorithm pDGM[) the estimates {zk,^^^^ \ka^^) satisfy p.4p . For this purpose we can 
choose the following values for fcout and Ein: 



(3.8) 



fout 



and Ein := 



1 



2(1 



Thus, we conclude from Theorems 13.21 13731 and 13^ that for these choices of fcout and 
£in the following estimates hold: 



/* - dp(Afc^„J < Eout, Zfeo„t e \\A'zk^^^ - b\\ < — £out and 
/3!|A*|| 9p \ (I ||Ao||2\ 

Finally, we are ready to summarize the above convergence rate analysis in the 
following algorithm. 

Algorithm [Inexact dual gradient method (IDGM)). 

Initialization: Choose parameters p > and < Ld < Choose an initial point 
Ao e K™. 

Outer iteration: For fc = 0, 1, . . . , fcout, perform: 

Step 1. {Inner loop). For given Afc, solve the inner problem p.ip with accu- 
racy £in, such that one of the stopping criterions ()2.4p - <{2.7\ are satisfied, to 
obtain Zk- 

Step 2. Form the approximate gradient vector of dp as Vdp{\k) '■— Azk — b. 
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Step 3. Select a^. G [i ,L^'^] and update := X]j=o '^r 
Step 4. Update Afe+i := \k + akVdp{Xk)- 
Output: := S^^^^ Ej=o "j^i- 

The penalty parameter p in this algorithm can be also updated adaptively by 
using e.g. the procedure given in [5]. 

3.2. Inexact dual fast gradient method. In this subsection we discuss a 



fast gradient scheme for solving the augmented Lagrangian dual problem (Dp). Fast 
gradient schemes were first proposed by Nesterov [15] and have also been discussed 
in the context of dual decomposition in [12]. A modification of these schemes for 
the case of inexact information can be also found in [3]. We shortly present such a 
scheme as follows. Given a positive sequence {Ok}k>o *- (0, +oo) with ~ Ij we 

define Sk '■= X]j=o assume that the sequence {dk}k>o satisfies Of^-^ = 5*^+1 

for all fc > 0. This condition leads to: 



(3.9) 0k+i := -(1 + J4el + 1) Vfc > and := 1. 



Note that the sequence {dk}k>o generated by (|3.9p is increasing and satisfies: 
(3.10) 0.5(fc + 1) < 6ifc < A; + 1 Vfc > 0. 

We can also obtain 0.25(fc + l)(fc + 2) < 5fc < 0.5(A: + l)(fc + 2) and Ej=o < 
(fc + l)(fc + 2)(fc + 3)/3. Now, we consider the dual fast gradient scheme as follows: 
given an initial point Ao G M™, we define two sequences of the dual variables {Xk}k>o 
and {/ifc}fc>o as: 



(IDFGM) 



Mfe -.^ Xk + L^^Vdp{Xk) 
Afe+i := (1 - ak+i)fik + flfc+i Ao + Y.i=o OiVdp{Xi) 



where the sequence ak+i '■= Sj^^^-^Ok+i- 

The following lemma, which represents an extension of the results in [12l [15] to 
the inexact case (see also [1]), will be used to derive estimates on both primal and 
dual suboptimality and also primal infeasibility for the proposed method. 

Lemma 3.5. \1S^ Under the assumptions of Theorem \ 3.1[ the two sequences 
{{Xk, t^k)}k>o generated by the dual fast gradient scheme (jlDFCMp satisfy: 

Skdp{iJik)> max J ['^"p(^^) + {^dp{X,), X - A,) ] - ^||A - Ao| 

- (l + y^Rp) Sin Vfc > 0. 

3=0 

The next theorem gives an estimate on dual suboptimality. 

Theorem 3.6. Under the assumptions of Theorem \3.1l let {{Xk, tJ'k)}k>o ^6 
the two sequences generated by the scheme (jIDFGM[) . Then, an estimate on dual 
suboptimality is given by the following expression: 
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Proof. By using inequality (|3.ip in p.lip wc obtain: 

Skdp{fJ.k) > Skdp{X*) - - AolP - CzSin^Sj. 

3=0 

Now, using the fact that Sk > 0.25(fc + l)(/s + 2) and ^^^q Sj < {k + l){k + 2){k + 3)/3 

and the definition of Cz = 1 + y^2LpRp, we obtain our resuh. □ 
We further define the following primal average sequence: 

k 

(3.12) 4:=5,7^^e,z,. 

j=o 

Next theorem gives an estimate on infeasibility of Zk for the original problem ((P|) . 

Theorem 3.7. Under the conditions of Theorem \3.6\ the point Zk defined by 
(|3.12p satisfies the following estimate on primal feasibility violation: 

(3.13) \\Azk-b\\<v[k,e,^), 

where v(k - ) - ^^^^^^ + / 2L,(k+3){i+^,R,)e.,^ 
wncie uyr.,^in) ■— (fc+i)(fc+2) ^ *\j 3(fc+l)(fe+2) 

Proof. By the definition of dp and Vc?p, the convexity of / and || ■ |p, and inequality 
(|3.12p we have: 

Y.e,[dp{\,) + (VJp(A,), A-A,)] ^J2^jfi-z,) + Sk {\,Azk-b)+Y,0,^\\Azj - bf 

j=0 j=0 j=0 

> Skfizk) + Sk {X,Azk - b)) + ^\\Azk - b\\\ 



dpilik) > /(^fc)+max {{X,Azk - b)) - ^o\\'} 



Substituting this inequality into (j3.1ip we obtain: 

k 

(3.14) + ^\\Azk - &f - CzSinS^' J2 

On the one hand, we can write: 

dp{^ik) - f{zk) - ^\\Azk ^bf< dp{\*) - f{zk) - ^\\Azk - bf 



(3.15) =min/:p(z,A*)-/(zfe)-|pzfe-6|p < {\\Azk~b). 



On the other hand, we have: 

(3.16) max| - ^j|A - Aof + (A, Az,. -b)\ = ^\\Azk - bf + (Aq, Az, - 6) 
Substituting ((XT5|) and (fXTC)) into ([5TTi)) we obtain: 

Q ^ 

-f\\Azk - bf + (Ao - \\Azk -b)< CzSinSk'Y.^^- 

^^•^ 3=0 
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If we define ^ := \\Azk — b\\, then the last mequahty hiiplies ^^-^-^^^^-^^'^ ^ Rd£, ^ 
4(fc+3) (g^g.^^ Therefore, we obtain ^ < i/(fc,£in), where •) is defined in p.l3p . □ 

Finally, we characterize the primal suboptimality for optimization problem ((P|) . 

Theorem 3.8. Under the conditions of Theorem \3.7\ the following estimates 
hold on primal suboptimality: 



Proof. The left-hand side inequality can be obtained similarly as in Theorem 
We now prove the right-hand side. From p.l4p and p.l6p we have: 

dpifik) > f{zk) + -bf + (Ao, Azk -b) + ^\\Azu - bf - Czs.nS-^ J2 

4(fc + 3) 

^^("^)' (fc + l)(fc + 2) "^"" 



Therefore, we get: 



2Ld n^ u2 . 4(fc-t-3)^ 



f^'^^ - '^^^^^ ^ (fc+i)(;+2) ""°" + 

Since dp[^k) < /*, we obtain the right-hand side inequality from the last relation. □ 
Similar to the previous subsection, we assume that we fix the outer accuracy Sout 
and the goal is to find fcout and a relation between Eont and ein such that after fcout 
outer iterations of the scheme pPFGMI) relations (|1.4I) holds. We can take e.g.: 



(3.17) 



8(i+v2i;:i?p)(fcout+3) 



Using now Theorems 13. 6[ 13.71 and 13.81 we can conclude that the following bounds for 
dual suboptimality, primal infeasibility and primal suboptimality hold: 

3 

/* - dp{\k,^^) < Eout, 2fco„t e Z, Pifco„t - ^11 < ^Eout and 

— 2^ 7 ' -/i^feoutj - J S ^ J ^out- 

We can summarize the above convergence rate analysis into the following algorithm. 
Algorithm {Inexact dual fast gradient method (IDFGM)). 

Initialization: Choose parameters p > and 6q := 1. Choose an initial point 

Ao e R™ and set 5*0 := 1. 

Outer iteration: For fc = 0, 1, . . . , fcout, perform: 

Step 1. {Inner loop). For given Afc, solve the inner problem p.ip with accu- 
racy £in, such that one of the stopping criterions ()2.4p - p.7p are satisfied, to 
obtain Zk. 

Step 2. Form the approximate gradient vector of dp as Vdp{Xk) := Azk — b. 
Step 3. Update pk ■= ^k + L^^Vdp{Xk). 

Step 4. Update O^+i := 0.5 (l + y/TTIelj , Sk+i := Sk + ^fe+i and Ok+i := 



fe+i- 
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Step 4. Update A^+i := (1 - ak+i)fJ.k + Ofc+i Ao + ^ Sj=o ^j Vdp(Aj 
Output: S^^l E'=o 

As in previous section, the penalty parameter p can be also updated adaptively 
by using the same procedure as before. 

4. Complexity certification for linear MPC problems. In this section we 
discuss different implementation aspects for the application of the algorithms derived 
in Sections 13.11 and 13 . 21 in the context of state-input constrained MPC for fast linear 
embedded systems. We first prove that for linear MPC with quadratic stage and final 
costs, the augmented Lagrangian function becomes strongly convex and therefore the 
inner problems (|2.1I) can be solved in linear time with a fast gradient scheme |14| . 
We also discuss how the different parameters, which appear in our derived complexity 
bounds of Algorithms pPGMp and pPFGMp . can be computed such that tight 
estimates on the total number of iterations can be derived and thus to facilitate the 
implementation on linear embedded systems with state-input constraints. 

4.1. Implementation aspects for MPC problems. Wc denote by Xjy a 
subset of the region of attraction for the MPC scheme discussed in Section 11.11 A 
detailed discussion on the stability of suboptimal MPC schemes can be found e.g. in 
[23] . For a given x € Xj\[, we denote with z*{x) an optimal solution for (P(a;)) and 
with A* (a;) an associated optimal multiplier. Usually, in MPC problems the stage and 
final costs are quadratic functions of the form: 

£{xi, Ui) := xfQxi + uj Rui and li{xN) '■— xJfPxN, 

where the matrices Q and P are positive semidefinite and R is positive definite. 
Note that in our formulation we do not require strongly convex stage cost, i.e. we 
do not impose the matrices Q and P to be positive definite. The following lemma 
characterizes the convexity properties of the augmented Lagrangian function. 



Lemma 4.1. // the optimization problem (P(a;) ) comes from a linear MPC prob- 
lem with quadratic stage and final costs, then the augmented Lagrangian Cp{z,X,x) is 
a strongly convex quadratic function w.r.t. variable z. 

Proof. If we consider quadratic costs in the MPC problem p.ip . then the objective 
function / is quadratic, i.e. f{z) := ^z'^Hz, where the Hessian H := diag((5,i?) is 
positive semidefinite, with Q := diag((3, • • • ,Q,P) and R := diag(i?, • • • ,R). Note 
that R is positive definite, since we assume R to be positive definite. Using these 
notations, we can rewrite the augmented Lagrangian in the form (see Section [TTT]) : 

Cp{z, A, x) := ^z'^iH + pA^A)z + {A^ \ ~ pA^h{x))z - b{xfX -f- ^h{xfh{x). 

It is straightforward to see that since H is positive semidefinite, then z"'" {H-\'pA^ A)z > 
for all z which satisfy Az ^ 0. On the other hand, if we consider the following set 
{z G M"| Az = 0}, which comes from the linear dynamics, we can rewrite equiva- 



lently this set as < z e M" | z = 



Au 
u 



entiC^}> wher 



, w e I 1-1 u >, wnere u := \uK ■ ■ ■ 



N-li 



and the matrix A is obtained from the matrices A^ and i?„ describing the dynamics 
of the system. Further, since Az = 0, we can write z^{H + pA^ A)z = z^Hz = 
A^QAu -\- u^Ru > for all u ^ 0. The last inequality follows from the fact that 



u 



R is positive definite. In conclusion, we proved that H + pA^ A is a positive definite 
matrix and therefore Cp{z, A, x) is a quadratic strongly convex function in 
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The previous lemma shows that in the hnear MPC case with quadratic costs the 
objective function of the inner subproblems £p are quadratic strongly convex in the 
first variable z. Moreover, Lp has also Lipschitz continuous gradient. Note that the 
convexity parameter cTp of this function can be computed easily: 

and the Lipschitz constant Lp of the gradient of Cp is: 

ip Amax(i^ + 

Note that since Cp{z^ A,.t) is strongly convex and with Lipschitz continuous gradient 
in the variable z, by solving the inner problem (|2.ip with a fast gradient scheme we can 
ensure stopping criterion (j2.4p in a linear number of inner iterations jl4j . Since the 
estimate for the number of inner iterations depends on (Tp, Lp and also on the diameter 
Rp of the set Z, we can see immediately that this diameter can be computed easily 
for cases when the set Z has a specific structure. Note that the set Z is a Cartesian 
product and thus: 



i?p ■.^^{N-l)Dl + Dl+NDl, 

where D^, D^^ and £>„ denotes the diameters of X, Xf and [/, respectively. These 
diameters can be computed explicitly for constraints sets defined e.g. by boxes or 
Euclidean balls, which typically appear in the context of MPC problems. 

Further, the estimates for the number of outer iterations depend on the norm of 
the dual optimal solution. We now discuss how we can bound ||A*|| in the MPC case. 
We make use of the result from [3] : 

Lemma 4.2. f^' For the family of MPC problems {P{x))x£Xn assume that 
there exists r > such that i?(0, r) C {Az — b{x) \ z € Z,x € X^^}, where i?(0, r) de- 
notes the Euclidean ball m 8^^^"^+""^ with center and radius r. Then, the following 
upper bound on the norm of the dual optimal solutions of MPC problems {P{x)) holds: 

Hi*. Nil ^ max^gz {Hz*{x),z- z*{x)) ^ ^ „ 

II A (x)|| < ; V.T e Xn, 

r 

where f := max {r | B{0, r) C {AZ - b{x) \ x € Xj^} } . 

Based on the previous lemma, in [TH] upper bounds are derived on ||A*(a;)|| for all 
X G Xpf for linear MPC problems with X, X{, U and Xn polyhedral sets: 

(4.1) 7^d > max ||A*(x)||. 

Recall that Lipschitz constant of the gradient of augmented dual function is Ld = 1/p. 

4.2. Total complexity of solving MPC problems. Now. wc assume that we 
know the outer accuracy £out and wc want to estimate the total number of iterations 
and also the number of flops per inner and outer iterations which have to be performed 
by Algorithms pPGMp or pPFGMp in order to solve the MPC problem (P(x)). 
For both algorithms we assume the initialization Ag = and the inner problems are 
solved using the stopping criterion (|2.4p . 



First, we discuss the complexity certificates in the case when problem (P(a;)) is 
solved using Algorithm pPGMp for all x £ X^. We denote by kf^ the number of 
inner iterations which has to be performed in order to solve each inner problem and 
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by fc^j the number of outer iterations. From the discussion in Section 13.11 an upper 
bound on the number of outer iterations is given by: 



(4.2) 



'''out ■ — 



^out 



Since we have proved that in the MFC case Cp{-,X,x) is strongly convex with con- 
vexity parameter cTp and has also Lipschitz continuous gradient with constant Lp, in 
order to find a point Zj,g (A) such that Cp{zf.G (A), A, x) — £p(z*(A), A, x) < ef^ we can 
apply a fast gradient scheme. From Theorem 2.2.3 in [T3] and taking into account 
that Ein = — „ V £out (see (|3.8p ) we get that the number of inner iterations for 



2(l + ^i?,p) 

finding such a point does not exceed 



(4.3) 




In the case of Algorithm pPFGMp . the number of outer iterations, according to 
the discussion in Section [3?2l is given by: 



(4.4) 



27^ri 



id 



Taking into account that in this case we consider that the inner accuracy is chosen 

Eout (see p.l7p ). then the number of inner iterations for 



as £ir] 



8(l+V^«p)Ku?+3) 

solving each inner problem will be given by 



(4.5) 



I.FG ._ 




2 / In ^ ^VLdT^d\/L^Rp (l 

f^n \ ^outv^^out 



Further, we are also interested in finding the total number of fiops for both outer 
and inner iterations. For solving the inner problem we use a simple fast gradient 
scheme for smooth strongly convex objective functions, see e.g. |14] . For this scheme, 
an inner iteration will require ti^°^'' :~ N (3n^ + 2nxnu + 2n\ + IQrix + 8n„) flops. 
Regarding the number of flops required by an outer iteration, the following values can 
be estabhshed: nl^'^ := N {2nl + 2n^,nu + Sn^,) + fc£nf„°P" for Algorithm (IIDGMI) 
and n^T'^^ — N {2nl + 2n^n„ + IQn^) + /Ci^'^nf"^' for Algorithm (lIDFGMp . re- 
spectively. 

5. Numerical experiments. In order to certify the efficiency of the proposed 
algorithms, we consider different numerical scenarios. We first analyze the behavior 
of Algorithms pPGMp and pPFGMp in terms of CPU time and number of itera- 
tions for some practical MFC problems and then we compare the CPU time, of our 
algorithms and of other well known QP solvers used in the context of MPC, on ran- 
domly generated QP problems. All the simulations were performed on a Laptop with 
CPU Intel T6670 with 2.2GHz and 4GB RAM memory, using Matlab R2008b. In aU 
simulations we consider Aq = 0. 

5.1. Practical MPC problems. In this section we apply the newly developed 
Algorithms pPGMp and pPFGMp on MPC problems for some practical applica- 
tions, i.e. a ball on plate system and an oscillating masses system. 
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5.1.1. Ball on plate system. The first application discussed in this section is 
the ball on plate system described in [T^. We consider box constraints for states X 
and Xf, inputs U and for the region of attraction Xm as in [12], while for the stage 
costs we take the matrices Q = qiqf , where qi = [2 1]'^ , R = 1 and we compute the 
terminal matrix P as the solution of the LQR problem. 

For different prediction horizons ranging from = 5 to A'^ = 20, we analyze first 
the behavior of Algorithms pDGM[) and pDFGM[) in terms of the number of outer 
iterations. For each prediction horizon length, we consider two different estimates for 
the number of outer iterations depending on the way we compute the upper bound on 
the optimal Lagrange multipliers \*{x). For Algorithm pPGMp . fc^j is the theoreti- 
cal number of iterations obtained using relation (j4.2[) with TZ^ computed according to 
[T^ (see (14.1^ ). while fc^^ ^^^^ is the average number of iterations obtained using our 
derived bound (j3.8p with i?d computed exactly using Gurobi 5.0.1 solver, iterations 
which correspond to 50 random initial states x £ Xm- We also compute the average 
number of outer iterations fc^^ ^.^^j observed in practice, obtained by imposing stop- 



ping criteria 



) - /* and \\Az^G^ 



b\\ less than accuracy Eout = 10 

FG „„;^„ ^ 

out real observed in practice. In all simulations we take p = 1. The 
results for both algorithms are reported in Figure lSTT] We can observe that in practice 



For Algorithm pPFGMp we compute in a similar way usmg 
using ((3T7)) and fcf <^ 



I.FG 
out, samp 



Algoritliiii IDGM 



Algoiitlim IDFCiM 



Prediction horizon N 



Fig. 5.1. Varia tion of kj ,, fc^,t,s.mp <^^<i ^Lt.real fi = {G;FG}) for Algorithm jlDGMj 
(left) and Algorithm l lIDFGMll (right) w.r.t the prediction horizon N, with accuracy Eout = 10~^. 

Algorithm pPFGMp performs much better than Algorithm pPGMp . Also, we can 
notice that the expected number of outer iterations /c^^. ^^^^^ and ^^^^ obtained 
from our derived bounds in Section [3] offer a good approximation for the real number 
of iterations performed by the two algorithms. Thus, these simulations show that 
our derived bounds in Section [3] are tight. On the other hand, is about three 
orders of magnitude, while fc^^. is about six orders of magnitude greater than the real 
number of iterations. 

In Figure [5?2] we also plot the evolution of the states and inputs over the simulation 
horizon for a prediction horizon = 5 and an outer accuracy £out = 10^'^. We 
observe that the system is driven to the equilibrium point. Since we obtained similar 
trajectories for the states and inputs with both algorithms, we present only the results 
for Algorithm pPFGMp . 

Since the number of outer iterations is also dependent on the way the inner 
accuracy £in is chosen, we are also interested in the behavior of the two algorithms 
w.r.t to Ein. For this purpose we apply Algorithms pPGMp and pPFGMp for 
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Fig. 5.2. The trajectories of the states and inputs for a prediction horizon N = 5 obtained 
using Algorithm l lIDFGMI I with accuracy Sout = lO""^. 



outer accuracy Eout = 10^"^ and varying e-m- In Figure [573] we plot the average number 
of outer iterations performed by the algorithms by taken 10 random samples for the 
initial state x € Xjy. We observe that we can increase the inner accuracy e^n derived in 




Inner accuracy "fn^ "in Inner accm-acy e^,. 



Fig. 5.3. The number of outer iterations performed by Algorithm llIDGMIl (left) and Algorithm 
llIDFGMIl (right) with fixed outer accuracy eout = 10^^ and different inner accuracies ein- 

Section |3] up to a certain value and the algorithms still perform a number of iterations 
less than the theoretical bounds derived in Section[3]for finding a suboptimal solution. 
On the other hand, if the inner accuracy is too large, the desired suboptimality cannot 
be ensured in a finite number of iterations. We see that Algorithm pPGMj) is less 
sensitive to the choice of the inner accuracy £in than Algorithm pPFGMp due to the 
fact that Algorithm pPFGMp accumulates errors (see Theorems 13.41 and 13. 8p . 

In conclusion, we notice from simulations that on the one hand the Algorithm 
pPFGMp is faster than pPGMp . but on the other hand that it is less robust. Thus, 
depending on the requirements of the application, one can choose between the two 
algorithms. 

5.1.2. Oscillating masses. The second example is a system comprised of M 
oscillating masses connected by springs to each other and to walls on either sides, 
having 2M states and M — 1 inputs. For a detailed description of the system, its 
parameters and constraints see [26] . We choose a quadratic stage cost with randomly 
generated positive semidefinite matrices Q E ]|j2A/x2A/^ having rank((5) = M, R = 
0.1 /m and the final cost P = Q. 

For this application we are interested in the CPU time. Thus, we consider only 
the Algorithm pPFGMp . which is usually faster than Algorithm pPGMp . In sim- 
ulations we vary the number M of masses and also the prediction horizon length A^. 
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Further, we consider both formulations of the MPC problem: sparse QP (i.e. we keep 
the states as variables) and dense QP (i.e. we eliminate the states using the dynamics 
of the system). Our goal is to compare the performances of Algorithm pDFGM|) 
and other methods used in the framework of linear MPC. Algorithm pPFGMp and 
Gurobi 5.0.1 solver (Gurl) arc used for solving the sparse formulation of the MPC 
problem. Alg. 1 in [17l[T6] and Gurobi 5.0.1 (Gur2) are used for solving the dense 
formulation of the MPC problem. In the implementation of the Algorithm (jIDFGMp 
we consider an adaptive scheme in order to update the penalty parameter p, similar 
to the one presented in [5| . Since the number of iterations is sensitive to the choice 
of penalty parameter, we have also tuned the initial guess of p. For each number of 
masses and prediction horizon, 50 simulations were run starting from different random 
initial states, We have considered the accuracy Cout = 10"'^ and the stopping criteria 
l/(%) ~ /*l cind \\Az]f — b\\ less than accuracy Eout- 



Table 5.1 

The average and maximum CPU time [s] (number of iterations) for Algorithm l lIDFGMII . Alg. 
1 in JjTy , Gurobi solver for sparse form (Gurl) and Gurobi solver for condensed form (Gur2). 



M 


N 


IDFGM 


Gurl 


Alg. 1 


Gur2 






avg 


max 


avg 


max 


avg 


max 


avg 


max 


5 


5 


0.03 (31) 


0.04 (33) 


0.007 (10) 


0.008 (11) 


0.05 (441) 


0.07 (604) 


0.005 (9) 


0.008 (11) 


5 


10 


0.07 (36) 


0.10 (51) 


0.009 (11) 


0.010 (12) 


0.13 (924) 


0.18 (1331) 


0.007 (12) 


0.008 (13) 


5 


20 


0.13 (65) 


0.22 (110) 


0.016 (11) 


0.017 (12) 


0.33 (1199) 


0.65 (2383) 


0.038 (12) 


0.043 (13) 


10 


5 


0.10 (28) 


0.12 (30) 


0.013 (10) 


0.014 (12) 


0.24 (1611) 


0.33 (2193) 


0.007 (10) 


0.008 (11) 


10 


10 


0.25 (47) 


0.38 (72) 


0.027 (11) 


0.028 (13) 


0.77 (2552) 


1.34 (4449) 


0.037 (12) 


0.041 (13) 


10 


20 


0.64 (70) 


1.25 (135) 


0.051 (12) 


0.055 (13) 


2.75 (2331) 


5.69 (4698) 


0.156 (10) 


0.168 (12) 


20 


5 


0.23 (42) 


0.34 (64) 


0.039 (11) 


0.043 (12) 


0.99 (3066) 


1.45 (4481) 


0.020 (11) 


0.025 (13) 


20 


10 


1.54 (98) 


2.90 (193) 


0.078 (12) 


0.084 (13) 


7.60 (5067) 


18.30 (11953) 


0.105 (12) 


0.115 (13) 


20 


20 


8.2 (356) 


14.9 (646) 


0.230 (12) 


0.770 (13) 


57.6 (12581) 


84.7 (18504) 


1.300 (12) 


2.120 (12) 



We can observe from Table [01 that Algorithm pDFGM[) outperforms Alg. 1 in 
[17j . especially when the dimension of the problem increases. On the other hand, we 
can notice that the solver Gurobi 5.0.1 performs much faster than our algorithm, 
since the MPC problem is sparse. However, the CPU times of the two algorithms are 
comparable in the case of dense QP problems (see next section). 

5.2. Random quadratic programming problems. In this section we com- 
pare the performance, in terms of CPU time, of Algorithms (jIDGMp and (jIDFGMp 
against some well known QP solvers used for solving MPC problems: quadprog (Mat- 
lab R2008b), Sedumi 1.3, Cplex 12.4 (IBM ILOG) and Gurobi 5.0.1. 

We consider random QP problems of the form 

min {0.5z'^Qz + q'^ z : s.t. ^2 = 6| , 

lb<z<ub 

where matrices Q € R""^" and A g RTtT^" are taken from a normal distribution 
with zero mean and unit variance. Matrix Q is then made positive semidefinite by 
transformation Q ^ Q^Q, having rank((5) ranging from 0.5n to 0.9n. Further, 
ub = —lb = 1 and b is taken from a uniform distribution. 

We plot in Figure [5^ the average CPU time for each solver, obtained by solving 
50 random QP's for each dimension n, with an accuracy Eout = 10"'^ and the stopping 
criteria \ f{zk) — f*\ and \\Azk — b\\ less than accuracy Eout- In the case of Algorithm 
(jIDGM[) . at each outer iteration we let the algorithm perform only k^-^ = 50 inner 
iterations. For the Algorithm pPFGMp we consider two scenarios: in the first 
one, we let the algorithm to perform only k^^'~^ = 100 inner iterations, while in the 
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second one we use the theoretic number of inner iterations obtained in Section [4.2l (see 
(|4.5|) ). As described previously, in our algorithms we consider an adaptive scheme for 
updating the penalty parameter p, similar to the one presented in [5]- We can observe 
that even if the Algorithms pPGMp and pPFGMp are well suited for embedded 
applications, i.e. the implementation of the iterates is very simple, the iteration 
complexity is low and also the number of iterations for finding an approximate solution 
can be easily predicted, the computational time is comparable with the one of the 
other solvers used in the context of MPC. Although the obtained averaged CPU times 
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Fig. 5.4. Average CPU time for solving QP problems of different sizes. 

are comparable, we cannot compare the exact computation complexity since for this 
purpose an equivalence between different stopping criteria of each solver should be 
studied, like e.g. the maximum violation of the constraints. 

6. Conclusions. Motivated by MPC problems for fast embedded linear sys- 
tems, we have proposed two dual gradient based methods for solving the augmented 
Lagrangian dual of a primal convex optimization problem with complicating linear 
constraints. We have moved the complicating constraints in the cost using augmented 
Lagrangian framework and solved the dual problem using gradient and fast gradient 
methods with inexact information. We have solved the inner subproblems only up to 
a certain accuracy, discussed the relations between the inner and the outer accuracy 
of the primal and dual problems and derived tight estimates on both primal and dual 
suboptimality and also on feasibility violation. We have also discussed some imple- 
mentation issues of the new algorithms for embedded linear MPC problems and tested 
them on several examples. 
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